% 2-D, MULTI-GROUP, k-eff IN CARTESIAN GEOMETRY
% j.roberts, 4/24/2010.
% I'm trying to clean up, verify, and speed up the diffusion solver

clear; 
%clc; %clf;

disp('-------------------------------------------')
disp('--     beginning 2D diffusion solver     --')
disp('-------------------------------------------')


%---- INPUT FILE (JUST AN "M" FILE) ---------------------------------------
twodbenchmark_input % 2-D IAEA benchmark, very close
%simple2d_input      % simple 1 box problem for reference
%simple2dRESP_input    % resp fct for 1 box problem
%twoelement_input % simple 2 box problem

tic

%---- GENERATE THE COEFFICIENTS -------------------------------------------
[AB, AL, AC, AR, AT, S, dx, dy, dv, dc, siga, sct] = ...
    twoDcoefMG( dat, abdat, numg, numm, xcm, xfm, ycm, yfm, mt, src, it, ....
    [BCL BCR BCB BCT], IBSL, IBSR, IBSB, IBST );

toc

nusig = S;
%---- SOLVE----------------------------------------------------------------
N = sum(xfm); M = sum(yfm);
% pad end of AB,AL
AB = [AB' zeros(numg,N+1)]';
AL = [AL' zeros(numg,1)]';
% pad top of AR,AT
AT = [zeros(numg,N+1) AT']';
AR = [zeros(numg,1) AR']';
%phi = zeros( (N+1)*(M+1), numg );

if it==0 %-----------------------------------------------------------------
%     do nothing
else % eigenvalue
    disp('begin keff iteration')
    nusig = S;  % just rename it right away to avoid confusion
    %s       = S;
    s = ones((N+1)*(M+1),1)/(ycm(end)+xcm(end));  % uniform guess
    keff    = 1.0;    % guess the initial keff
    errK    = 1;      errS    = 1; 
    epsK    = 1e-9;   epsS = 1e-3;
    iter    = 0;      itmx = 300;
    phi     = zeros((N+1)*(M+1),numg); 
    A1 = spdiags([AB(:,1) AL(:,1) AC(:,1) AR(:,1) AT(:,1)], ...
                    [-N-1 -1 0 1 N+1],(N+1)*(M+1),(N+1)*(M+1));
    A2 = spdiags([AB(:,2) AL(:,2) AC(:,2) AR(:,2) AT(:,2)], ...
                    [-N-1 -1 0 1 N+1],(N+1)*(M+1),(N+1)*(M+1));

    while  ( ( errK>epsK ) || (errS>epsS) ) && iter < itmx
        % set to solve 1st group flux
        s = s/sum(s);
        phi(:,1) = A1 \ (s(:,1)/keff) ;  % assumes all fiss n's fast
        if numg > 1
            for g = 2:numg
                scsrc   = zeros( (N+1)*(M+1), 1);
                % compute scattering source
                for gg = 1:g-1 
                    scsrc(:,1) = scsrc(:,1) + sct(:,g,gg).*phi(:,gg);
                end
                phi(:,g) = A2 \ ( scsrc(:,1) ); % rhs is scatter only
            end   
        end
        sold = s; 
        kold = keff;   
        s(:,1) = 0; % reset fission source
        for g = 1:numg
            s(:,1) = s(:,1) + nusig(:,g).*phi(:,g);
        end         
        keff = sum(s.*dv)*kold/sum(sold.*dv);
        indx = find(s);
        errS = max( abs((s(indx)-sold(indx))./s(indx)) );
        errK = abs( (keff-kold)/keff );
        iter = iter + 1;  % number of iterations
        if ( mod(iter,20)== 0)
            disp([' iter = ',num2str(iter), ' keff = ',num2str(keff)])
            disp([' errK = ',num2str(errK), ' errS = ',num2str(errS)])
        end
    end
    disp([ 'final result: keff = ',num2str(keff)])
    disp([ '              iter = ',num2str(iter)])
    keff
end

toc

max(phi(:,1)./phi(:,2))


